Analysis of Sudan survey data 2012 - 2013

This R-markdown document contains all r-code used to carry out the analysis for the paper on fish catch rates and species densities (species richness) along the Red Sea coast of Sudan. The analysis is based on three surveys: November 2012, May 2013 and November 2013.

Here the analysis is structured around the 7 management regions of the Sudanese Red Sea coast as follows:

All code and data are stored on GitHub

Libraries, data etc.

Libraries, dependencies and functions

Load and manipulate data

Reading catch, station and traits data.

Neither station data, nor catch data has a complete depth record, but by combining depth para from each we can get a complete depth parameter for all stations.

Sudan map and management areas

Loads the map data for Sudan, loads the management areas from shape-file etc.

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/eriko/Documents/GitHub/Sudan2019/Sudan-master/sudan_management_areas", layer: "sudan_regions"
## with 7 features
## It has 3 fields
## Integer64 fields read as strings:  id

Allocating catch positions to management areas

Adds ‘Area’ to each line in the catch table.

Modifying dataset and creating one dataset for CPUE analysis and another for traits

Of the catch data, 12 fish registrations lack weight, (i.e. this was forgotten entered into the database during the survey). These registrations, would if included lead to 12 more registrations of CPUE = 0.

Adding traits to ‘catch’ data set (without 0-catch stations)

Some more data wrangling that adds the traits from the traits table to the catch data, using only the station with catches (there are no traits for ‘NOCATCH’ species).

Also, only select traps, gillnets and handlines as these were the only gear with sufficient numbers and consistent use to be analyzed.

Lastly adds number of gear deployed at each station to each line.

Bathymetric map of Sudan

Fig. 1 in MS

A nice, bathymetric map of Sudan with management areas overlaid.

Currently fonts and colours are adjusted for use on the ICES ASC 2019 poster, removing the bathy-contour lines (se ‘lwd’ to 0), and varying the colours of management area text to white for areas 1-4.

## quartz_off_screen 
##                 2

Stations plotted on maps, pr survey

Faceted map plotting position of all catch stations for each of the three surveys.

## Saving 7 x 5 in image

Station information

Station table

Table describing the sampling effort in each area pr survey, number of different gear types, max, min and average depth, number of traps with / without catches.

survey id Ntraps Nhl NGn TBhrs HLhrs GNhrs DepthAvg DepthMax DepthMin
2012901 1 22 0 0 694.066 0.000 0.000 42.45455 142 13
2012901 2 54 3 3 721.549 78.000 62.000 40.89744 71 0
2012901 3 26 0 1 451.183 0.000 14.000 30.80769 95 8
2012901 4 5 0 0 77.084 0.000 0.000 22.60000 54 10
2012901 5 31 0 4 677.896 0.000 78.317 21.32258 30 7
2012901 6 36 0 1 850.417 0.000 38.250 32.38889 88 0
2012901 7 31 0 8 712.200 0.000 135.849 31.06667 66 5
2013002 1 29 0 1 420.163 0.000 24.000 31.48148 70 5
2013002 2 81 3 5 1102.379 377.500 221.569 27.09091 145 5
2013002 3 32 0 1 545.811 0.000 24.000 28.87500 60 0
2013002 4 13 0 10 160.266 0.000 137.815 29.50000 67 9
2013002 5 33 2 5 208.899 192.000 195.766 20.46154 50 7
2013002 6 45 1 2 642.413 156.000 78.000 34.70270 88 9
2013002 7 39 2 0 666.410 186.000 0.000 33.46154 76 5
2013005 1 23 1 4 271.551 3.000 84.000 29.82353 80 10
2013005 2 57 2 6 500.101 111.000 156.000 38.27273 80 7
2013005 3 9 0 2 123.099 0.000 36.000 26.50000 70 9
2013005 4 16 2 4 151.184 4.500 142.767 32.90000 68 12
2013005 5 30 2 3 317.498 7.000 146.947 25.52381 65 11
2013005 6 40 4 2 443.983 31.501 71.915 40.14815 89 6
2013005 7 22 0 2 171.784 0.000 203.171 34.75000 54 11

Species table

Number of fish (organized by family and species) caught by Gillnet or Traps for each survey, and in total across all surveys and gears.

  • (new table for revision 2020) *
    fam_name Sci_name 2012901_GN 2012901_TB 2013002_GN 2013002_TB 2013005_GN 2013005_TB sum
    ACANTHURIDAE Acanthurus nigricans 0 10 0 24 6 49 89
    ACANTHURIDAE Acanthurus nigricauda 0 0 0 12 0 0 12
    SOLEIDAE Achirus sp. * 0 0 0 0 1 0 1
    HOLOCENTRIDAE Adioryx ruber 0 0 0 11 0 4 15
    HOLOCENTRIDAE Adioryx spinifer 0 0 0 0 1 0 1
    SERRANIDAE Aethaloperca rogaa 0 2 0 10 0 1 13
    ALBULIDAE Albula vulpes 0 0 0 0 5 0 5
    CARANGIDAE Alectis indicus 0 0 0 0 1 0 1
    CARANGIDAE Alepes vari 0 0 0 0 4 0 4
    SPARIDAE Argyrops filamentosus 0 0 0 7 0 0 7
    SPARIDAE Argyrops sp. 0 20 0 0 0 0 20
    SPARIDAE Argyrops spinifer 0 0 0 6 0 5 11
    ARIIDAE Arius heudelotii 0 6 0 0 0 0 6
    ARIIDAE Arius thalassinus 0 0 0 2 1 0 3
    SCOMBRIDAE Auxis thazard 14 0 0 0 0 0 14
    BALISTIDAE Balistapus undulatus 0 0 0 0 0 1 1
    BALISTIDAE Balistoides viridescens 0 0 0 1 0 0 1
    BOTHIDAE Bothus pantherinus 0 0 0 0 1 0 1
    CAESIONIDAE Caesio caerulaurea 0 0 15 0 0 0 15
    CAESIONIDAE Caesio suevicus 0 0 0 0 1 0 1
    CARANGIDAE Carangoides armatus 0 0 0 0 4 0 4
    CARANGIDAE Carangoides bajad 13 1 6 14 49 0 83
    CARANGIDAE Carangoides ferdau 0 0 0 1 7 0 8
    CARANGIDAE Carangoides fulvoguttatus 0 0 0 0 16 0 16
    CARANGIDAE Carangoides sp. 0 0 0 0 1 0 1
    CARANGIDAE Caranx (Caranx) melampygus 2 0 0 1 6 0 9
    CARANGIDAE Caranx (Caranx) sexfasciatus 16 0 22 2 43 0 83
    CARANGIDAE Caranx (Gnathanodon) speciosus 0 0 0 0 4 0 4
    CARANGIDAE Caranx ignobilis 0 0 0 1 1 0 2
    CARANGIDAE Caranx sp. 0 0 0 0 1 0 1
    CARCHARHINIDAE Carcharhinus albimarginatus 2 0 0 0 0 0 2
    S H A R K S Carcharhinus melanopterus 9 0 0 0 5 0 14
    S H A R K S Carcharhinus wheeleri 0 0 1 0 1 0 2
    SERRANIDAE Cephalopholis argus 0 1 0 0 0 0 1
    SERRANIDAE Cephalopholis miniatus ? 0 0 0 0 1 0 1
    CHAETODONTIDAE Chaetodon auriga 0 3 0 0 0 0 3
    CHAETODONTIDAE Chaetodon semilarvatus 0 0 16 2 0 0 18
    CHANIDAE Chanos chanos 0 0 0 0 1 0 1
    LABRIDAE Cheilinus lunulatus 0 0 0 1 0 0 1
    LABRIDAE Cheilinus quinquecintus 0 1 0 0 0 0 1
    CHIROCENTRIDAE Chirocentrus dorab 15 0 17 0 60 0 92
    CARANGIDAE Chorinemus lysan 0 0 0 6 0 0 6
    PLATYCEPHALIDAE Cociella crocodila 0 0 0 0 1 0 1
    MUGILIDAE Crenimugil crenilabis 0 0 1 0 0 0 1
    CARANGIDAE Decapterus macarellus 0 0 1 0 0 0 1
    SCOMBRIDAE Decapterus russelli 0 0 1 4 0 0 5
    HAEMULIDAE Diagramma pictum 0 1 0 0 0 0 1
    DIODONTIDAE Diodon hystrix 0 0 0 0 3 0 3
    ECHENEIDIDAE Echeneis naucrates 0 0 0 3 2 0 5
    SCOMBRIDAE Elagatis bipinnulata 0 0 6 0 0 0 6
    SERRANIDAE Epinephelus chlorostigma 0 0 0 1 0 0 1
    SERRANIDAE Epinephelus fuscoguttatus 0 7 0 7 1 2 17
    SERRANIDAE Epinephelus sexfasciatus 0 0 0 1 0 1 2
    SERRANIDAE Epinephelus summana 0 0 2 1 0 0 3
    SERRANIDAE Epinephelus tauvina 1 11 4 4 3 2 25
    SCOMBRIDAE Euthynnus affinis 0 0 4 0 3 0 7
    JUVENILES FISHJUV 9 0 0 1 0 0 10
    GERREIDAE Gerres oyena 0 0 1 0 6 0 7
    SCOMBRIDAE Grammatorcynus bicarinatus 0 0 3 0 0 0 3
    SCOMBRIDAE Grammatorcynus bilineatus 26 2 2 0 5 0 35
    LETHRINIDAE Gymnocranius microdon 0 0 1 0 0 0 1
    LETHRINIDAE Gymnocranius robinsoni 0 0 0 2 3 0 5
    SCOMBRIDAE Gymnosarda unicolor 0 0 0 0 10 0 10
    MURAENIDAE Gymnothorax flavimarginatus 0 3 0 0 0 0 3
    MURAENIDAE Gymnothorax javanicus 0 28 0 11 0 1 40
    HAEMULIDAE Hectromynicus pictus 0 0 0 4 0 0 4
    HEMIRAMPHIDAE Hemiramphus far 0 0 0 0 1 0 1
    SCARIDAE Hipposcarus harid 0 0 4 0 4 0 8
    SCOMBRIDAE Katsuwonus pelamis 1 0 0 0 0 0 1
    KYPHOSIDAE Kyphosus vaigiensis 0 0 0 0 10 0 10
    LETHRINIDAE Lethrinus elongatus * 0 10 1 17 7 5 40
    LETHRINIDAE Lethrinus harak 0 0 0 0 10 0 10
    LETHRINIDAE Lethrinus lentjan 3 33 6 47 31 28 148
    LETHRINIDAE Lethrinus mahsena 0 16 2 19 0 18 55
    LETHRINIDAE Lethrinus mahsenoides * 0 4 0 0 6 0 10
    LETHRINIDAE Lethrinus nebulosus 0 0 0 0 0 1 1
    LETHRINIDAE Lethrinus ramak 0 0 1 1 0 1 3
    LETHRINIDAE Lethrinus xanthochilus 0 0 0 2 0 1 3
    LUTJANIDAE Lutjanus argentimaculatus 0 0 0 1 0 1 2
    LUTJANIDAE Lutjanus bohar 0 32 10 69 3 8 122
    LUTJANIDAE Lutjanus cf fluviflamma 0 0 0 0 2 0 2
    LUTJANIDAE Lutjanus coccineus * 0 0 0 0 0 1 1
    LUTJANIDAE Lutjanus ehrenbergii 6 0 11 0 7 0 24
    LUTJANIDAE Lutjanus fulviflamma 0 0 0 0 2 0 2
    LUTJANIDAE Lutjanus gibbus 0 38 0 51 1 28 118
    LUTJANIDAE Lutjanus kasmira 0 4 0 3 0 2 9
    LUTJANIDAE Lutjanus monostigma 2 4 0 3 0 1 10
    LUTJANIDAE Lutjanus rivulatus 0 0 0 1 0 0 1
    LUTJANIDAE Lutjanus sebae 0 0 0 1 0 0 1
    LUTJANIDAE Lutjanus sp. 0 1 0 0 0 0 1
    LUTJANIDAE Macolor niger 0 0 0 2 1 0 3
    LETHRINIDAE Monotaxis grandoculis 0 0 1 0 0 0 1
    MULLIDAE Mulloides flavolineatus 0 0 1 0 0 0 1
    MULLIDAE Mulloides vanicolensis 0 0 0 0 1 0 1
    HOLOCENTRIDAE Myripristis murdjan 0 0 3 0 4 0 7
    ACANTHURIDAE Naso hexacanthus 0 0 18 18 37 0 73
    ACANTHURIDAE Naso lituratus 0 0 0 0 1 0 1
    NO CATCH NO CATCH 3 0 0 0 0 0 3
    LUTJANIDAE Paracaesio sordius 0 2 0 0 0 0 2
    EPHIPPIDAE Platax boersi 0 1 0 0 0 0 1
    EPHIPPIDAE Platax orbicularis 0 2 0 2 1 2 7
    HAEMULIDAE Plectorhinchus gaterinus 0 7 1 0 3 0 11
    POMADASYIDAE (HAEMULIDAE) Plectorhinchus pictus 0 0 0 0 1 0 1
    HAEMULIDAE Plectorhynchus pictus 0 0 0 0 1 0 1
    HAEMULIDAE Plectorhynchus schotaf 0 0 0 0 1 0 1
    SERRANIDAE Plectropomus pessuliferus 0 2 0 2 0 0 4
    PRIACANTHIDAE Priacanthus hamrur 0 0 4 0 0 0 4
    LUTJANIDAE Pristipomoides multidens 0 9 0 7 0 0 16
    BALISTIDAE Pseudobalistes flavimarginatus 0 1 0 0 0 0 1
    SCOMBRIDAE Rastrelliger kanagurta 0 0 6 0 21 0 27
    SCOMBRIDAE Sarda orientalis 0 0 1 0 0 0 1
    HOLOCENTRIDAE Sargocentron spiniferum 0 13 1 13 0 5 32
    SCARIDAE Scarus ferrugineus 0 0 1 0 0 0 1
    SCARIDAE Scarus frenatus 0 0 2 0 2 0 4
    SCARIDAE Scarus ghobban 0 0 0 0 1 0 1
    SCOMBRIDAE Scomber japonicus 0 0 0 1 0 0 1
    CARANGIDAE Scomberoides lysan 15 20 60 18 45 0 158
    CARANGIDAE Scomberoides tol 0 0 4 0 74 0 78
    SCOMBRIDAE Scomberomorus commerson 35 0 0 0 0 0 35
    SCOMBRIDAE Scomberomorus lineolatus 0 0 0 0 9 0 9
    SCOMBRIDAE Scomberomorus tritor 2 0 0 0 0 0 2
    SIGANIDAE Siganus argenteus 0 0 3 0 0 0 3
    SIGANIDAE Siganus luridus 1 0 1 0 0 0 2
    SIGANIDAE Siganus rivulatus 0 0 0 0 1 0 1
    SIGANIDAE Siganus stellatus 0 0 1 0 0 2 3
    SPARIDAE Sparus sp. 0 0 0 6 0 0 6
    SPHYRAENIDAE Sphyraena forsteri 0 0 1 0 0 0 1
    SPHYRAENIDAE Sphyraena jello 3 0 0 0 0 0 3
    SPHYRAENIDAE Sphyraena putnamie 0 0 0 0 1 0 1
    SPHYRAENIDAE Sphyraena qenie 0 0 6 0 4 0 10
    CARCHARHINIDAE Sphyrna lewini 1 0 0 0 0 0 1
    R A Y S Taeniura lymma 0 0 0 0 1 0 1
    SCOMBRIDAE Thunnus albacares 11 0 0 0 0 0 11
    CARCHARHINIDAE Triaenodon obesus 0 3 0 10 0 0 13
    BELONIDAE Tylosurus choram 2 0 4 0 0 0 6
    MUGILIDAE Valamugil engeli 0 0 0 0 1 0 1
    SERRANIDAE Variola louti 0 0 0 2 0 0 2

Number of set traps per depth

Violin plot of depths at wich traps were set. I prefer this visualization of depth ranges of the traps, rather than a (boring) table.

## quartz_off_screen 
##                 2

Test differneces in depths in each area between years and between areas

## 
##  Shapiro-Wilk normality test
## 
## data:  c3$depth
## W = 0.90427, p-value < 2.2e-16

## 
##  Kruskal-Wallis rank sum test
## 
## data:  c3$depth and c3$id
## Kruskal-Wallis chi-squared = 36.877, df = 6, p-value = 1.861e-06
## Warning in posthoc.kruskal.nemenyi.test.default(c3$depth, c3$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  c3$depth and c3$id 
## 
##   1      2       3      4      5      6     
## 2 1.0000 -       -      -      -      -     
## 3 0.8599 0.6487  -      -      -      -     
## 4 0.9804 0.9477  1.0000 -      -      -     
## 5 0.0025 1.2e-05 0.3022 0.4513 -      -     
## 6 0.9943 0.9561  0.9883 0.9997 0.0075 -     
## 7 0.9983 0.9860  0.9835 0.9994 0.0117 1.0000
## 
## P value adjustment method: none
## [1] "2012901"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 29.371, df = 6, p-value = 5.174e-05
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1      2      3      4      5      6     
## 2 0.9991 -      -      -      -      -     
## 3 0.9060 0.4375 -      -      -      -     
## 4 0.6436 0.3884 0.9600 -      -      -     
## 5 0.1017 0.0011 0.7700 1.0000 -      -     
## 6 0.7221 0.1350 1.0000 0.9811 0.8655 -     
## 7 0.9290 0.4572 1.0000 0.9415 0.6429 0.9995
## 
## P value adjustment method: none 
## [1] "2013002"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 13.422, df = 6, p-value = 0.0368
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1     2     3     4     5     6    
## 2 0.999 -     -     -     -     -    
## 3 0.997 1.000 -     -     -     -    
## 4 1.000 1.000 1.000 -     -     -    
## 5 0.405 0.460 0.783 0.768 -     -    
## 6 0.999 0.893 0.918 0.999 0.078 -    
## 7 1.000 0.989 0.988 1.000 0.226 1.000
## 
## P value adjustment method: none 
## [1] "2013005"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$depth and ci$id
## Kruskal-Wallis chi-squared = 10.394, df = 6, p-value = 0.109
## Warning in posthoc.kruskal.nemenyi.test.default(ci$depth, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$depth and ci$id 
## 
##   1    2    3    4    5    6   
## 2 1.00 -    -    -    -    -   
## 3 1.00 0.98 -    -    -    -   
## 4 1.00 1.00 1.00 -    -    -   
## 5 0.64 0.13 0.99 0.85 -    -   
## 6 1.00 0.98 1.00 1.00 0.66 -   
## 7 1.00 1.00 1.00 1.00 0.68 1.00
## 
## P value adjustment method: none

Depth data were not normally distributed and hence the difference in depth could only be analyzed using non-parametric Kruskal - Wallis rank sum test.

Of the 63 survey - area comparisons, depth distributions were different between areas as follows: Nov.2012 survey: 1:5, 2:5, 5:6 and 5:7
May 2013 survey: 2:5

No significant differences for the Nov. 2013 survey

GAM model of trap depths

Fits a GAM model (using GAMLSS package) to the depth with area as explanatory variable and survey as a random factor. Use ‘chooseDIst’ function to test which distribution on the positive real line best fit the data.

## GAMLSS-RS iteration 1: Global Deviance = 5840.466 
## GAMLSS-RS iteration 2: Global Deviance = 5840.466
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## minimum GAIC(k= 2 ) family: exGAUS 
## minimum GAIC(k= 3.84 ) family: exGAUS 
## minimum GAIC(k= 6.51 ) family: exGAUS
##                 2     3.84     6.51
## EXP            NA       NA       NA
## GA             NA       NA       NA
## IG             NA       NA       NA
## LOGNO          NA       NA       NA
## LOGNO2         NA       NA       NA
## WEI            NA       NA       NA
## WEI2           NA       NA       NA
## WEI3           NA       NA       NA
## IGAMMA         NA       NA       NA
## PARETO2  6076.280 6091.000 6112.360
## PARETO2o 6028.005 6042.725 6064.085
## GP             NA       NA       NA
## BCCG           NA       NA       NA
## BCCGo          NA       NA       NA
## exGAUS   5676.647 5693.207 5717.237
## GG             NA       NA       NA
## GIG            NA       NA       NA
## LNO            NA       NA       NA
## BCTo           NA       NA       NA
## BCT            NA       NA       NA
## BCPEo          NA       NA       NA
## BCPE           NA       NA       NA
## GB2            NA       NA       NA
## GAIG with k= 6.51
##   exGAUS PARETO2o  PARETO2      EXP       GA       IG    LOGNO   LOGNO2 
## 5717.237 6064.085 6112.360       NA       NA       NA       NA       NA 
##      WEI     WEI2     WEI3   IGAMMA       GP     BCCG    BCCGo       GG 
##       NA       NA       NA       NA       NA       NA       NA       NA 
##      GIG      LNO     BCTo      BCT    BCPEo     BCPE      GB2 
##       NA       NA       NA       NA       NA       NA       NA
## GAIG with k= 6.51 
## GAMLSS-RS iteration 1: Global Deviance = 5702.051 
## GAMLSS-RS iteration 2: Global Deviance = 5682.573 
## GAMLSS-RS iteration 3: Global Deviance = 5669.27 
## GAMLSS-RS iteration 4: Global Deviance = 5662.855 
## GAMLSS-RS iteration 5: Global Deviance = 5660.216 
## GAMLSS-RS iteration 6: Global Deviance = 5659.202 
## GAMLSS-RS iteration 7: Global Deviance = 5658.841 
## GAMLSS-RS iteration 8: Global Deviance = 5658.714 
## GAMLSS-RS iteration 9: Global Deviance = 5658.675 
## GAMLSS-RS iteration 10: Global Deviance = 5658.656 
## GAMLSS-RS iteration 11: Global Deviance = 5658.65 
## GAMLSS-RS iteration 12: Global Deviance = 5658.648 
## GAMLSS-RS iteration 13: Global Deviance = 5658.647
## GAIG with k= 6.51 
## ******************************************************************
## Family:  c("exGAUS", "ex-Gaussian") 
## 
## Call:  gamlss(formula = depth ~ factor(id) + random(factor(survey)),  
##     family = names(getOrder(t1, 3)[1]), data = c3) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15.15720    1.47773  10.257   <2e-16 ***
## factor(id)2  0.08425    1.66309   0.051   0.9596    
## factor(id)3 -3.55885    2.00296  -1.777   0.0761 .  
## factor(id)4 -4.87739    2.41340  -2.021   0.0437 *  
## factor(id)5 -3.79141    1.82704  -2.075   0.0384 *  
## factor(id)6 -3.68338    1.77533  -2.075   0.0384 *  
## factor(id)7 -2.10235    1.89522  -1.109   0.2677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7951     0.0917   19.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  log 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.92573    0.05178    56.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  9.000015
##       Residual Deg. of Freedom:  665 
##                       at cycle:  13 
##  
## Global Deviance:     5658.647 
##             AIC:     5676.647 
##             SBC:     5717.266 
## ******************************************************************

## 
## Family:  c("exGAUS", "ex-Gaussian") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = depth, family = "exGAUS", data = c3) 
## 
## Mu Coefficients:
## [1]  12.55
## Sigma Coefficients:
## [1]  1.753
## Nu Coefficients:
## [1]  2.955
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   671 
## Global Deviance:     5673.47 
##             AIC:     5679.47 
##             SBC:     5693.01

The GAMLSSS model picked ‘exGAUS’ as the best distribution to use in the model, which seems sensible from plotting the model over historgram of the depth.

Using the ‘exGAUS’ model areas 3, 4, 5 and 6 come out as significant factors, meaning the depth in these areas are significantly different from depths in other areas. However, the model does not say which pairwise area comparisons (within each survey) are significantly different.

Combination of GAMLSS model and Kruskal - Wallis

The GAMLSS model shows that there are signicant effects of area 3, 4, 5, and 6 on the depth distribution, and the K-W test shows significant differences for the Nov. 2012 survey between area 1:5, 2:5, 5:6, 5:7 and for the May 2013 survey for areas 2:4, which is in corrspondence with the GAM results.

## GAMLSS-RS iteration 1: Global Deviance = 5839.315 
## GAMLSS-RS iteration 2: Global Deviance = 5839.315
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## minimum GAIC(k= 2 ) family: exGAUS 
## minimum GAIC(k= 3.84 ) family: exGAUS 
## minimum GAIC(k= 6.51 ) family: exGAUS
## GAIG with k= 6.51
##   exGAUS PARETO2o  PARETO2      EXP       GA       IG    LOGNO   LOGNO2 
## 5727.751 6076.768 6125.020       NA       NA       NA       NA       NA 
##      WEI     WEI2     WEI3   IGAMMA       GP     BCCG    BCCGo       GG 
##       NA       NA       NA       NA       NA       NA       NA       NA 
##      GIG      LNO     BCTo      BCT    BCPEo     BCPE      GB2 
##       NA       NA       NA       NA       NA       NA       NA
## GAIG with k= 6.51 
## GAMLSS-RS iteration 1: Global Deviance = 5700.529 
## GAMLSS-RS iteration 2: Global Deviance = 5680.767 
## GAMLSS-RS iteration 3: Global Deviance = 5667.163 
## GAMLSS-RS iteration 4: Global Deviance = 5660.614 
## GAMLSS-RS iteration 5: Global Deviance = 5657.787 
## GAMLSS-RS iteration 6: Global Deviance = 5656.739 
## GAMLSS-RS iteration 7: Global Deviance = 5656.404 
## GAMLSS-RS iteration 8: Global Deviance = 5656.231 
## GAMLSS-RS iteration 9: Global Deviance = 5656.178 
## GAMLSS-RS iteration 10: Global Deviance = 5656.152 
## GAMLSS-RS iteration 11: Global Deviance = 5656.145 
## GAMLSS-RS iteration 12: Global Deviance = 5656.142 
## GAMLSS-RS iteration 13: Global Deviance = 5656.141 
## GAMLSS-RS iteration 14: Global Deviance = 5656.141
## GAIG with k= 6.51 
## ******************************************************************
## Family:  c("exGAUS", "ex-Gaussian") 
## 
## Call:  gamlss(formula = depth ~ factor(id) + factor(survey),  
##     family = names(getOrder(t1, 3)[1]), data = c3) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            15.0137     1.6905   8.881   <2e-16 ***
## factor(id)2             0.1809     1.6538   0.109   0.9129    
## factor(id)3            -3.3130     1.9969  -1.659   0.0976 .  
## factor(id)4            -4.9122     2.3832  -2.061   0.0397 *  
## factor(id)5            -3.7182     1.8163  -2.047   0.0410 *  
## factor(id)6            -3.7938     1.7696  -2.144   0.0324 *  
## factor(id)7            -2.0126     1.8879  -1.066   0.2868    
## factor(survey)2013002  -0.6828     1.1398  -0.599   0.5493    
## factor(survey)2013005   1.0894     1.2113   0.899   0.3688    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.77954    0.09392   18.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  log 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.92830    0.05183    56.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  11
##       Residual Deg. of Freedom:  663 
##                       at cycle:  14 
##  
## Global Deviance:     5656.141 
##             AIC:     5678.141 
##             SBC:     5727.786 
## ******************************************************************

Running the GAMLSS model with survey as a factor (not random), shows that the November 2012 survey has significant effects on the depth, while the other two surveys don’t have a significant effect on depth. The conclusion is therefore that depth distribution of traps is significantly different in for the Nov. 2012 survey compared to the other two, with areas 3,4,5, and 6 being significantly different on the area level.

Duration of trap sets

Mean, min & max hrs. of fishing

id maxhrs minhrs meanhrs
1 46.417 11.500 20.48641
2 44.500 7.917 16.47669
3 24.250 8.167 16.86333
4 27.033 13.033 16.11076
5 40.200 11.667 17.79179
6 30.083 8.000 19.00436
7 31.866 10.417 19.47891
##        ss                  id             fhrs       
##  Min.   :2.013e+09   Min.   :1.000   Min.   : 7.917  
##  1st Qu.:2.013e+09   1st Qu.:2.000   1st Qu.:15.425  
##  Median :2.013e+09   Median :4.000   Median :16.150  
##  Mean   :2.013e+09   Mean   :3.909   Mean   :17.984  
##  3rd Qu.:2.013e+09   3rd Qu.:6.000   3rd Qu.:18.700  
##  Max.   :2.013e+09   Max.   :7.000   Max.   :46.417

Fishing hours (traps) by survey and area

Boxplot of median with mean values plotted as red circle. The average trap soak-time (fishing hours) across all surveys and areas was 17.984 hrs with a median of 16.15, but the November 2012 survey had more varied and higher soak times than the May 2013 and Nov. 2013.

Test differences in soak times (traps) K-W test

## 
##  Shapiro-Wilk normality test
## 
## data:  c3$Fhrs
## W = 0.67297, p-value < 2.2e-16

## 
##  Kruskal-Wallis rank sum test
## 
## data:  c3$Fhrs and c3$id
## Kruskal-Wallis chi-squared = 68.971, df = 6, p-value = 6.645e-13
## Warning in posthoc.kruskal.nemenyi.test.default(c3$Fhrs, c3$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  c3$Fhrs and c3$id 
## 
##   1       2       3       4       5       6      
## 2 0.89835 -       -       -       -       -      
## 3 1.00000 0.88439 -       -       -       -      
## 4 0.61831 0.95548 0.60026 -       -       -      
## 5 0.99997 0.65429 1.00000 0.42582 -       -      
## 6 0.14066 1.3e-05 0.20277 0.00167 0.19214 -      
## 7 0.01625 2.0e-07 0.02891 0.00012 0.02234 0.97141
## 
## P value adjustment method: none
## [1] "2012901"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 60.806, df = 6, p-value = 3.087e-11
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1       2       3       4       5       6      
## 2 0.07376 -       -       -       -       -      
## 3 0.06785 0.99940 -       -       -       -      
## 4 0.18427 0.94372 0.98714 -       -       -      
## 5 0.96585 0.48177 0.39946 0.46876 -       -      
## 6 0.99285 0.00023 0.00085 0.04610 0.49836 -      
## 7 0.96104 8.6e-05 0.00033 0.02854 0.32747 0.99990
## 
## P value adjustment method: none 
## [1] "2013002"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 42.845, df = 6, p-value = 1.252e-07
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1       2       3       4       5       6      
## 2 0.99993 -       -       -       -       -      
## 3 0.23422 0.14820 -       -       -       -      
## 4 0.99997 1.00000 0.74873 -       -       -      
## 5 0.37056 0.28089 0.99998 0.85136 -       -      
## 6 0.00475 0.00024 0.94518 0.20902 0.83861 -      
## 7 0.06190 0.01735 0.99966 0.51167 0.99399 0.99505
## 
## P value adjustment method: none 
## [1] "2013005"
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ci$Fhrs and ci$id
## Kruskal-Wallis chi-squared = 27.538, df = 6, p-value = 0.0001148
## Warning in posthoc.kruskal.nemenyi.test.default(ci$Fhrs, ci$id, "Chisq"): Ties
## are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  ci$Fhrs and ci$id 
## 
##   1     2     3     4     5     6    
## 2 0.375 -     -     -     -     -    
## 3 0.307 0.966 -     -     -     -    
## 4 0.962 0.995 0.870 -     -     -    
## 5 0.302 1.000 0.995 0.973 -     -    
## 6 0.979 0.830 0.627 1.000 0.719 -    
## 7 0.958 0.015 0.045 0.496 0.015 0.426
## 
## P value adjustment method: none

Of the 63 area- survey combinations 16 were significantly different (Kruskal Wallis rank sum test with Nemenyi post-hoc test, p<0.05). In Nov. 2012 area 1 was significantly different from 7, area 2 from 6 and7, area 3:7, area 4 vs 6 and 7, area 5 vs 6. In May 2013 area 2 was significantly different from 6 and 7, area 3 from 6 and 7, and area 4 from 6 and 7. In Nov 2013 area 1 was significantly different from area 6, while area 2 was significantly different from area 6 and 7.

GAMLSS model of soak times

## GAMLSS-RS iteration 1: Global Deviance = 3903.587 
## GAMLSS-RS iteration 2: Global Deviance = 3903.587
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations

## Warning in additive.fit(x = X, y = wv, w = wt * w, s = s, who = who,
## smooth.frame, : additive.fit convergence not obtained in 30 iterations
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## minimum GAIC(k= 2 ) family: BCPEo 
## minimum GAIC(k= 3.84 ) family: BCPEo 
## minimum GAIC(k= 6.51 ) family: BCPEo
## GAIG with k= 6.51
##    BCPEo     BCPE     BCTo      BCT      GB2       IG    BCCGo   exGAUS 
## 3257.667 3265.534 3304.446 3314.150 3407.192 3536.386 3559.944 3566.905 
##     BCCG       GG   IGAMMA      GIG    LOGNO      LNO       GA   LOGNO2 
## 3567.756 3570.111 3573.235 3579.743 3597.631 3597.631 3630.957 3686.328 
##      WEI     WEI3     WEI2      EXP PARETO2o       GP  PARETO2 
## 3788.868 3788.868 3803.669 5291.946 5331.031 5360.799 5367.490
## GAIG with k= 6.51 
## GAMLSS-RS iteration 1: Global Deviance = 3522.762 
## GAMLSS-RS iteration 2: Global Deviance = 3209.979 
## GAMLSS-RS iteration 3: Global Deviance = 3189.191 
## GAMLSS-RS iteration 4: Global Deviance = 3182.858 
## GAMLSS-RS iteration 5: Global Deviance = 3181.535 
## GAMLSS-RS iteration 6: Global Deviance = 3180.58 
## GAMLSS-RS iteration 7: Global Deviance = 3179.169 
## GAMLSS-RS iteration 8: Global Deviance = 3178.542 
## GAMLSS-RS iteration 9: Global Deviance = 3176.934 
## GAMLSS-RS iteration 10: Global Deviance = 3175.871 
## GAMLSS-RS iteration 11: Global Deviance = 3175.791 
## GAMLSS-RS iteration 12: Global Deviance = 3174.872 
## GAMLSS-RS iteration 13: Global Deviance = 3174.879 
## GAMLSS-RS iteration 14: Global Deviance = 3173.936 
## GAMLSS-RS iteration 15: Global Deviance = 3173.325 
## GAMLSS-RS iteration 16: Global Deviance = 3173.287 
## GAMLSS-RS iteration 17: Global Deviance = 3173.038 
## GAMLSS-RS iteration 18: Global Deviance = 3173.037
## GAIG with k= 6.51 
## ******************************************************************
## Family:  c("BCPEo", "Box-Cox Power Exponential-orig.") 
## 
## Call:  gamlss(formula = Fhrs ~ factor(id) + random(factor(survey)),  
##     family = names(getOrder(t1, 3)[1]), data = c3) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  2.834797   0.001474 1922.615  < 2e-16 ***
## factor(id)2 -0.014547   0.001879   -7.740 3.74e-14 ***
## factor(id)3 -0.018907   0.003349   -5.646 2.45e-08 ***
## factor(id)4 -0.034830   0.005583   -6.239 7.88e-10 ***
## factor(id)5 -0.006294   0.003503   -1.797   0.0729 .  
## factor(id)6  0.070165   0.003430   20.458  < 2e-16 ***
## factor(id)7  0.032570   0.001915   17.012  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.2200     0.1469  -8.303 5.73e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.5773     0.4479  -5.754 1.33e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.81374    0.08432  -9.651   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  12.99999
##       Residual Deg. of Freedom:  661 
##                       at cycle:  18 
##  
## Global Deviance:     3173.037 
##             AIC:     3199.037 
##             SBC:     3257.709 
## ******************************************************************
## Warning in MLE(ll4, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma, :
## possible convergence problem: optim gave code=1 false convergence (8)

## 
## Family:  c("BCPEo", "Box-Cox Power Exponential-orig.") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = Fhrs, family = "BCPEo", data = c3) 
## 
## Mu Coefficients:
## [1]  2.767
## Sigma Coefficients:
## [1]  -1.54
## Nu Coefficients:
## [1]  -1.495
## Tau Coefficients:
## [1]  -0.4741
## 
##  Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   670 
## Global Deviance:     3448.48 
##             AIC:     3456.48 
##             SBC:     3474.53

GAMLSS found the ‘BCPEo’ to give the best fit to the ‘Fhrs’ data (see histogram plot).

Using the ‘BCPEo’ distribution a model of Fhrs by area (as factor) and survey as a random factor showed that all areas had a significant effect on the soak time for traps, meaning that soak time varied siginificantly between the areas.

If including survey as a factor in the model all surveys also had a significant effect of Fhrs, showing that soak time varied significantly between all surveys and areas.

Analysis of catches by traits and other factors

Catch rates pr area & survey

Must aggregate data per station

## quartz_off_screen 
##                 2

Statistical analysis CPUE catches

Testing for normality of CPUEw of trap catches

Shapiro test for normality and Q-Q plots shows that data is non-normal and zero-inflated.

## 
##  Shapiro-Wilk normality test
## 
## data:  cpue.tr.st$CPUEw
## W = 0.55028, p-value < 2.2e-16
## 
## Attaching package: 'qqplotr'
## The following objects are masked from 'package:ggplot2':
## 
##     stat_qq_line, StatQqLine

### Test for difference in CPUEw of traps between areas (non-parametric test) Both CPUEw is not significantly different between the areas,

Significant difference for CPUEW between Nov 2013 survey and the Nov 2012 and May 2013 surveys, but not between the Nov12 and May13 surveys

## 
##  Kruskal-Wallis rank sum test
## 
## data:  cpue.tr.st$CPUEw and cpue.tr.st$id
## Kruskal-Wallis chi-squared = 10.307, df = 6, p-value = 0.1123
## Warning in posthoc.kruskal.nemenyi.test.default(cpue.tr.st$CPUEw,
## cpue.tr.st$id, : Ties are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  cpue.tr.st$CPUEw and cpue.tr.st$id 
## 
##   1    2    3    4    5    6   
## 2 1.00 -    -    -    -    -   
## 3 1.00 1.00 -    -    -    -   
## 4 0.99 0.99 1.00 -    -    -   
## 5 0.94 0.90 0.98 1.00 -    -   
## 6 0.95 0.79 0.89 0.70 0.22 -   
## 7 1.00 1.00 1.00 1.00 1.00 0.55
## 
## P value adjustment method: none
## 
##  Kruskal-Wallis rank sum test
## 
## data:  cpue.tr.st$CPUEw and cpue.tr.st$survey
## Kruskal-Wallis chi-squared = 33.528, df = 2, p-value = 5.241e-08
## Warning in posthoc.kruskal.nemenyi.test.default(cpue.tr.st$CPUEw,
## cpue.tr.st$survey, : Ties are present. Chi-sq was corrected for ties.
## 
##  Pairwise comparisons using Nemenyi-test with Chi-squared    
##                        approximation for independent samples 
## 
## data:  cpue.tr.st$CPUEw and cpue.tr.st$survey 
## 
##         2012901 2013002
## 2013002 0.98    -      
## 2013005 9.1e-06 7.5e-07
## 
## P value adjustment method: none

Zero-inflated model (better approach)

This method is more appropriate than the Kruskal-Wallis approach above as it evaluates the effects of all factors at the same time, instead of evaluating the effects of area separate from survey.

First, developed a model for CPUE (weight and numbers) including depth, area, survey and gear (only traps and gillnets), which performed better than a simpler model excluding gear.

Also tested models where survey was included as a random factor (using two different approaches: ´random´and ´re´), although this excludes evaluation of the potential significant effects of the surveys on the results.

Tested models for both CPUE by weight (CPUEw) and CPUE by numbers (CPUEn).

Model performance were compared based on AIC.

## GAMLSS-RS iteration 1: Global Deviance = 882.3705 
## GAMLSS-RS iteration 2: Global Deviance = 882.3705
## GAMLSS-RS iteration 1: Global Deviance = 883.5556 
## GAMLSS-RS iteration 2: Global Deviance = 883.5553
## GAMLSS-RS iteration 1: Global Deviance = 883.5553 
## GAMLSS-RS iteration 2: Global Deviance = 883.555
## GAMLSS-RS iteration 1: Global Deviance = 824.5923 
## GAMLSS-RS iteration 2: Global Deviance = 824.5951 
## GAMLSS-RS iteration 3: Global Deviance = 824.595
## GAMLSS-RS iteration 1: Global Deviance = 891.7485 
## GAMLSS-RS iteration 2: Global Deviance = 891.7485
##            df      AIC
## mod4 23.29925 871.1935
## mod  13.00000 908.3705
## mod3 12.73369 909.0224
## mod2 12.73357 909.0224
## mod5 13.72109 919.1906
## GAMLSS-RS iteration 1: Global Deviance = 2954.59 
## GAMLSS-RS iteration 2: Global Deviance = 2954.59
## GAMLSS-RS iteration 1: Global Deviance = 2955.38 
## GAMLSS-RS iteration 2: Global Deviance = 2955.379
## GAMLSS-RS iteration 1: Global Deviance = 622.7372 
## GAMLSS-RS iteration 2: Global Deviance = 622.7371
## GAMLSS-RS iteration 1: Global Deviance = 2653.48 
## GAMLSS-RS iteration 2: Global Deviance = 2653.48
##             df       AIC
## modA3 13.00000  648.7371
## mod4  23.29925  871.1935
## mod   13.00000  908.3705
## mod3  12.73369  909.0224
## mod2  12.73357  909.0224
## modA4 14.00000 2681.4797
## modA1 14.00000 2982.5900
## modA2 14.20733 2983.7937
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEw ~ depth + factor(id) + factor(gear) +  
##     random(factor(survey)), nu.formula = ~depth + factor(id) +  
##     factor(gear) + random(factor(survey)), family = ZAGA,      data = cpue.st) 
## 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.196781   0.217732   0.904    0.366    
## depth          -0.002022   0.002785  -0.726    0.468    
## factor(id)2     0.219134   0.191534   1.144    0.253    
## factor(id)3     0.027183   0.240497   0.113    0.910    
## factor(id)4    -0.359623   0.269631  -1.334    0.183    
## factor(id)5    -0.256229   0.226000  -1.134    0.257    
## factor(id)6    -0.162825   0.204251  -0.797    0.426    
## factor(id)7    -0.191843   0.220544  -0.870    0.385    
## factor(gear)TB -1.488538   0.172201  -8.644   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.02777    0.03243   0.856    0.392
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.612486   0.399613  -4.035 6.05e-05 ***
## depth          -0.003167   0.004074  -0.777    0.437    
## factor(id)2     0.182223   0.270936   0.673    0.501    
## factor(id)3     0.255953   0.334275   0.766    0.444    
## factor(id)4     0.385587   0.394158   0.978    0.328    
## factor(id)5     0.453593   0.312951   1.449    0.148    
## factor(id)6    -0.412901   0.297750  -1.387    0.166    
## factor(id)7     0.307303   0.308661   0.996    0.320    
## factor(gear)TB  1.730137   0.335438   5.158 3.24e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  734 
## Degrees of Freedom for the fit:  23.29925
##       Residual Deg. of Freedom:  710.7007 
##                       at cycle:  3 
##  
## Global Deviance:     824.595 
##             AIC:     871.1935 
##             SBC:     978.3353 
## ******************************************************************

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.01075036 
##                        variance   =  1.031094 
##                coef. of skewness  =  0.008704951 
##                coef. of kurtosis  =  3.290931 
## Filliben correlation coefficient  =  0.9986442 
## ******************************************************************

## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEn ~ depth + factor(id) + factor(survey) +  
##     factor(gear), family = ZAGA, data = cpue.st) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -0.111560   0.212432  -0.525   0.5996    
## depth                 -0.005006   0.002471  -2.026   0.0432 *  
## factor(id)2            0.343525   0.167665   2.049   0.0408 *  
## factor(id)3            0.076865   0.213978   0.359   0.7195    
## factor(id)4            0.190577   0.234250   0.814   0.4162    
## factor(id)5            0.467716   0.198977   2.351   0.0190 *  
## factor(id)6           -0.216463   0.178400  -1.213   0.2254    
## factor(id)7            0.246084   0.191877   1.283   0.2001    
## factor(survey)2013002  0.205185   0.109024   1.882   0.0602 .  
## factor(survey)2013005  0.284145   0.132690   2.141   0.0326 *  
## factor(gear)TB        -1.840150   0.151532 -12.144   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.14492    0.03367  -4.304 1.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.04360    0.07384   0.591    0.555
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  734 
## Degrees of Freedom for the fit:  13
##       Residual Deg. of Freedom:  721 
##                       at cycle:  2 
##  
## Global Deviance:     622.7371 
##             AIC:     648.7371 
##             SBC:     708.5178 
## ******************************************************************

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.01702982 
##                        variance   =  0.979929 
##                coef. of skewness  =  0.433784 
##                coef. of kurtosis  =  4.98197 
## Filliben correlation coefficient  =  0.9868307 
## ******************************************************************

Model for weight of catch (CPUEw)

Gear was the only factor significantly affecting CPUEw in any of the models run, also in model 4 (with the lowest AIC).

Depth had no significant effect on CPUE in these models, neither did area.

Model for number of catch (CPUEn)

Depth, area (2 & 5), surveys and gear were all found to have a significant effect on CPUEn.

The model with survey as a normal factor (not random) had the lowest AIC. This also allowed to evaluate wether surveys had significant effects on the CPUE.

GAM model CPUE - traps only

Statistical analysis of plot of trap catch rates

## GAMLSS-RS iteration 1: Global Deviance = 680.1101 
## GAMLSS-RS iteration 2: Global Deviance = 680.1094
## GAMLSS-RS iteration 1: Global Deviance = 679.8366 
## GAMLSS-RS iteration 2: Global Deviance = 679.8366
## GAMLSS-RS iteration 1: Global Deviance = 395.6543 
## GAMLSS-RS iteration 2: Global Deviance = 395.6522 
## GAMLSS-RS iteration 3: Global Deviance = 395.6523
## GAMLSS-RS iteration 1: Global Deviance = 394.265 
## GAMLSS-RS iteration 2: Global Deviance = 394.2647
##             df      AIC
## mod55 12.00000 418.2647
## mod44 11.77455 419.2014
## mod33 12.00000 703.8366
## mod22 12.66004 705.4294
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEn ~ depth + factor(id) + factor(survey),  
##     family = ZAGA, data = c22) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -2.1132812  0.1940001 -10.893  < 2e-16 ***
## depth                 -0.0008637  0.0026047  -0.332  0.74030    
## factor(id)2            0.2826635  0.1679995   1.683  0.09294 .  
## factor(id)3            0.1774132  0.2149544   0.825  0.40947    
## factor(id)4            0.0003821  0.2591314   0.001  0.99882    
## factor(id)5            0.5122892  0.2040252   2.511  0.01228 *  
## factor(id)6           -0.2127678  0.1785497  -1.192  0.23383    
## factor(id)7            0.1951116  0.1968309   0.991  0.32192    
## factor(survey)2013002  0.2953217  0.1088671   2.713  0.00685 ** 
## factor(survey)2013005  0.2910664  0.1377350   2.113  0.03496 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.20337    0.03668  -5.544 4.29e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.16154    0.07752   2.084   0.0376 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  670 
## Degrees of Freedom for the fit:  12
##       Residual Deg. of Freedom:  658 
##                       at cycle:  2 
##  
## Global Deviance:     394.2647 
##             AIC:     418.2647 
##             SBC:     472.3521 
## ******************************************************************

Using CPUE-weight data for traps in a simlar ZAGA GAMLSSS model (with survey as a random factor) showed that neither area or depth had significant effects on the CPUEw of traps. Including survey as an ordinary factor in the model (mod33) increased the fit (lower AIC), and identified that the Nov 2013 survey had significant effect on the CPUEw - that the observed lower catch rates during this survey were significantly different from the other two surveys.

Similar models using CPUE-numbers had lower AIC, with the model using survey as an ordinary factor (mod55) having the lowest AIC. For mod55 all surveys, as well as areas 2 and 5 were significant factors.

## 
##  Welch Two Sample t-test
## 
## data:  subset(cpue.st, gear == "TB")$depth and (subset(cpue.st, gear == "GN")$depth)
## t = 3.5683, df = 67.859, p-value = 0.0006648
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   6.034582 21.347788
## sample estimates:
## mean of x mean of y 
##  31.73806  18.04688

Depth was, however, not a significant factor in any of the four models developed for the trap CPUE, by weight or numbers. The depth distribution of bottom depth at the trap and gillnet stations (Gillnet fishing depth was at the surface, while trap fishing depth = bottom depth) were significantly different (t-test, p<0.05, df=67.9 - see also frequncy plot above). Therefore, the models including both gear types led to a wider span and a more diverse distribution of the depths than for the traps-only data set.

Combine plots of traps in one panels

Use ‘patchwork’ package to arrange plots

## Saving 7 x 5 in image

CPUE by gear type - by family and trophic group

Only for gillnets and traps

## Warning: Ignoring unknown parameters: binwidt

## Warning: Ignoring unknown parameters: binwidt

Test differences in CPUE By Species in trap catches, by range of traits and factors

Use GAM model and evaluate how CPUE by weight (CPUEw) was affected by physical / area /surve factors and the traits : - depth - survey (random factor) - area - family group - trophic group - trophic level (linked to trophic group) - Place in water column - Diel activity - Habitat - Gregariousness - MaxLength

## GAMLSS-RS iteration 1: Global Deviance = -1214.863 
## GAMLSS-RS iteration 2: Global Deviance = -1214.762 
## GAMLSS-RS iteration 3: Global Deviance = -1214.761
## GAMLSS-RS iteration 1: Global Deviance = -1215.724 
## GAMLSS-RS iteration 2: Global Deviance = -1215.724
## GAMLSS-RS iteration 1: Global Deviance = 966.4375 
## GAMLSS-RS iteration 2: Global Deviance = 966.4373
## GAMLSS-RS iteration 1: Global Deviance = 966.4372 
## GAMLSS-RS iteration 2: Global Deviance = 966.4368
##             df       AIC
## mod8  25.00000 -1165.724
## mod7  25.51665 -1163.728
## mod9  25.00019  1016.438
## mod10 26.00000  1018.437
## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEn ~ depth + factor(id) + factor(survey) +  
##     TrophicLevel + factor(TGShort) + factor(WaterCol) +  
##     factor(DielActivity) + factor(Habitat) + factor(Gregariousness) +  
##     MaxLength, family = ZAGA, data = na.omit(subset(c4,      gear == "TB"))) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                0.260653   0.769440   0.339 0.734978
## depth                                      0.002477   0.001700   1.457 0.145878
## factor(id)2                               -0.115302   0.112872  -1.022 0.307646
## factor(id)3                               -0.032661   0.139215  -0.235 0.814637
## factor(id)4                               -0.111354   0.239163  -0.466 0.641768
## factor(id)5                                0.043209   0.135927   0.318 0.750744
## factor(id)6                               -0.344605   0.121391  -2.839 0.004770
## factor(id)7                               -0.060450   0.128857  -0.469 0.639246
## factor(survey)2013002                      0.225175   0.063262   3.559 0.000418
## TrophicLevel                               0.059661   0.151419   0.394 0.693792
## factor(TGShort)Coral.                     -0.177659   0.425717  -0.417 0.676680
## factor(TGShort)Herb.                       0.459963   0.322578   1.426 0.154713
## factor(TGShort)Invert.                    -0.167304   0.106783  -1.567 0.117995
## factor(TGShort)Plankt.                     3.451168   0.695630   4.961 1.06e-06
## factor(WaterCol)Demersal                  -2.364709   0.373519  -6.331 6.83e-10
## factor(WaterCol)pelagic non-site attached -0.773423   0.309798  -2.497 0.012961
## factor(WaterCol)pelagic site attached     -4.920564   0.914149  -5.383 1.28e-07
## factor(DielActivity)Night                 -0.302945   0.129815  -2.334 0.020130
## factor(Habitat)Coral                       2.318828   1.007965   2.301 0.021956
## factor(Habitat)sand                        0.116507   0.579656   0.201 0.840811
## factor(Gregariousness)2                   -0.027771   0.073892  -0.376 0.707247
## factor(Gregariousness)3                   -1.031894   0.352940  -2.924 0.003664
## MaxLength                                 -0.009479   0.001407  -6.737 5.96e-11
##                                              
## (Intercept)                                  
## depth                                        
## factor(id)2                                  
## factor(id)3                                  
## factor(id)4                                  
## factor(id)5                                  
## factor(id)6                               ** 
## factor(id)7                                  
## factor(survey)2013002                     ***
## TrophicLevel                                 
## factor(TGShort)Coral.                        
## factor(TGShort)Herb.                         
## factor(TGShort)Invert.                       
## factor(TGShort)Plankt.                    ***
## factor(WaterCol)Demersal                  ***
## factor(WaterCol)pelagic non-site attached *  
## factor(WaterCol)pelagic site attached     ***
## factor(DielActivity)Night                 *  
## factor(Habitat)Coral                      *  
## factor(Habitat)sand                          
## factor(Gregariousness)2                      
## factor(Gregariousness)3                   ** 
## MaxLength                                 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.56730    0.03408  -16.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.1355     0.2478  -12.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  408 
## Degrees of Freedom for the fit:  25
##       Residual Deg. of Freedom:  383 
##                       at cycle:  2 
##  
## Global Deviance:     -1215.724 
##             AIC:     -1165.724 
##             SBC:     -1065.442 
## ******************************************************************

… this needs more thinking….

Evaluting catches of fish in relation to their traits. Should be evaluated by numbers, for single large fish not to dominate the results. The model based on CPUE by numbers (mod8) had much lower AIC than models with numbers and fishing hours as a factor. Also, including survey as an ordinary factor led to lower AIC than including it as a random factor.

Model (mod8) results show that numbers of fish caught (CPUE by numbers) was significantly affected by: - area 6 - Nov. 2013 survey - trophic level - trophic group: invertivores and planktivores - place in water column (all 3 classes) - diel activity (night) - coral habitat - gregariousness: 2&3 - Max length

Interpretation of results: There is no effect of depth, wich is consistent with the models for CPUEn above (on all data, without the traits). The Nov 2013 survey and area 6 were significant factors, indicating that these were differnt from the other areas & surveys. Surprisingly, the trophic groups invertivores and planktivores were significant factors. Place in water colum was significant for all categories of this trait, nocturnal fish prefering coral habitats, siz of the fish, and of gregariouslevels 2 & 3 were also significant factors.

Biological explanation (possible): These results indicates that catch in numbers of traps is dependent on demersal and other site-attached fish, linked to coral habitats, with some gregariousness, of certain max-length in an area. This is consistent with traps being set at the bottom close to coral reefs (not on sandy habitats), and the way traps work - attracting gregarious fish with bait. The entry into the traps is limited to fish of a certain size range, and the mesh size of the traps means that only fish larger than a certain size were retained.

The significance of invertivores and planktivores in relation to the catches is more difficult to interpret, but these are often species occurring in schools attracted to debris of bait spread out by larger fish feeding on the trap bait.

Also, area 6 and the November 2013 survey were significant factors pointing to systematic differnences in the catch rates for this area and survey, consistent with the results from the CPEUn analysis of all gear stations (without including traits) above.

CPUE by trophic group by area & by survey (year)

FIG 5 in Feb 2019 version of MS

Presents the CPUE pr hr, for each management area and survey, standardized by the number of traps set in each area X survey combination.

Trap catches were dominated by carnivores, occuring in all area X survey combinations, followed by invertivores who also occurred everywhere, albeight in lower densities than the carnivores.

Cathces of planktivores and invertivores varied substantially between areas and surveys, in contrast to carnivores who were caught in all surveys and areas. This explains why invertivores and planktivores were signifant factors on the catches (GAM model of CPUEn with traits), while carnivores were not.

Area 2, 5 and 6 seem to have the highest diversity of trophic groups, but this has to be calculated as functional diversity, using the FD package (or similar).

## quartz_off_screen 
##                 2

Trophic group CPUE by depth

A modified version of Figure 6 in Feb 2019 MS.

Basically a more complex version of the previous plot. Adds the depth dimension to the plot, and the surveys as colours.

Shows uniform catch rates of carnivores and invertivores from 0-50m.

All surveys seem to overlap, so no difference in trophich group depth dependent catch rates between the surveys.

## quartz_off_screen 
##                 2

Trophic level in Trap catches

## quartz_off_screen 
##                 2

Maxlength in catches

## quartz_off_screen 
##                 2

Place in water column

## quartz_off_screen 
##                 2

Diel activity

## quartz_off_screen 
##                 2

Gregariousness

## quartz_off_screen 
##                 2

Biodiversity (species densities, trophic diversity etc)

Aggregated Species density per area for (traps only)

Species per area pr hour of fishing, traps only.

Area Nov.2012 May 2013 Nov.2013 Average
1 0.0094519 0.0196548 0.0214763 0.0147296
2 0.0121593 0.0134669 0.0121096 0.0127123
3 0.0193318 0.0252437 0.0246737 0.0231029
4 0.0389186 0.0295256 0.0250687 0.0286485
5 0.0140647 0.0185659 0.0132908 0.0152575
6 0.0143005 0.0138886 0.0142919 0.0141270
7 0.0148488 0.0138863 0.0171116 0.0148840

The highest number of species was observed in area 2 (36), followed by area 3, 6 and 7 (25 in each).

Average species density across the three surveys(number of species pr hour of fishing - soak time of traps) ranged from 0.01271 (area 2) to 0.028 (area 4). Area 4 had the highest species density across all surveys (table above), while the area with the lowest species density varied by survey: area 1 in November 2012 , area 2 in May 2013 , and area 3 in November 2013.

Species density pr station - Zero Adjusted GAM model

## GAMLSS-RS iteration 1: Global Deviance = -4285.522 
## GAMLSS-RS iteration 2: Global Deviance = -4285.522
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function
## Warning in body(fun): argument is not a function
## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged

## Warning in RS(): Algorithm RS has not yet converged
## minimum GAIC(k= 2 ) family: BCPE 
## minimum GAIC(k= 3.84 ) family: BCPE 
## minimum GAIC(k= 6.51 ) family: BCPE
##                  2      3.84      6.51
## EXP      -2463.854 -2445.454 -2418.754
## GA       -4163.753 -4143.513 -4114.143
## IG       -4047.663 -4027.423 -3998.053
## LOGNO    -4135.944 -4115.704 -4086.334
## LOGNO2          NA        NA        NA
## WEI      -4085.168 -4064.928 -4035.558
## WEI2     -4084.146 -4063.906 -4034.536
## WEI3     -4085.168 -4064.928 -4035.558
## IGAMMA   -4099.195 -4078.955 -4049.585
## PARETO2   6450.000  6470.240  6499.610
## PARETO2o -2388.665 -2368.425 -2339.055
## GP       -2399.313 -2379.073 -2349.703
## BCCG     -4196.277 -4174.197 -4142.157
## BCCGo    -4187.901 -4165.821 -4133.781
## exGAUS   -4349.439 -4327.359 -4295.319
## GG       -4175.515 -4153.435 -4121.395
## GIG      -4161.753 -4139.673 -4107.633
## LNO      -4135.944 -4115.704 -4086.334
## BCTo     -4452.069 -4428.149 -4393.439
## BCT      -4463.401 -4439.481 -4404.771
## BCPEo    -4507.118 -4483.198 -4448.488
## BCPE     -4519.402 -4495.482 -4460.772
## GB2      -4358.274 -4334.354 -4299.644
## GAIG with k= 6.51
##      BCPE     BCPEo       BCT      BCTo       GB2    exGAUS      BCCG     BCCGo 
## -4460.772 -4448.488 -4404.771 -4393.439 -4299.644 -4295.319 -4142.157 -4133.781 
##        GG        GA       GIG     LOGNO       LNO    IGAMMA       WEI      WEI3 
## -4121.395 -4114.143 -4107.633 -4086.334 -4086.334 -4049.585 -4035.558 -4035.558 
##      WEI2        IG       EXP        GP  PARETO2o   PARETO2    LOGNO2 
## -4034.536 -3998.053 -2418.754 -2349.703 -2339.055  6499.610        NA
## GAIG with k= 6.51 
## GAMLSS-RS iteration 1: Global Deviance = -4423.958 
## GAMLSS-RS iteration 2: Global Deviance = -4514.165 
## GAMLSS-RS iteration 3: Global Deviance = -4528.121 
## GAMLSS-RS iteration 4: Global Deviance = -4530.868 
## GAMLSS-RS iteration 5: Global Deviance = -4534.342 
## GAMLSS-RS iteration 6: Global Deviance = -4537.106 
## GAMLSS-RS iteration 7: Global Deviance = -4537.984 
## GAMLSS-RS iteration 8: Global Deviance = -4540.997 
## GAMLSS-RS iteration 9: Global Deviance = -4541.756 
## GAMLSS-RS iteration 10: Global Deviance = -4543.27 
## GAMLSS-RS iteration 11: Global Deviance = -4543.821 
## GAMLSS-RS iteration 12: Global Deviance = -4544.342 
## GAMLSS-RS iteration 13: Global Deviance = -4544.655 
## GAMLSS-RS iteration 14: Global Deviance = -4544.881 
## GAMLSS-RS iteration 15: Global Deviance = -4545.111 
## GAMLSS-RS iteration 16: Global Deviance = -4545.187 
## GAMLSS-RS iteration 17: Global Deviance = -4545.319 
## GAMLSS-RS iteration 18: Global Deviance = -4545.353 
## GAMLSS-RS iteration 19: Global Deviance = -4545.404 
## GAMLSS-RS iteration 20: Global Deviance = -4545.402
## Warning in RS(): Algorithm RS has not yet converged
## GAIG with k= 6.51
## Warning in summary.gamlss(fm): summary: vcov has failed, option qr is used instead
## ******************************************************************
## Family:  c("BCPE", "Box-Cox Power Exponential") 
## 
## Call:  gamlss(formula = Species_Hrs ~ depth + factor(area) + factor(survey),  
##     family = names(getOrder(t1, 3)[1]), data = spue.df) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                         Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)            5.108e-02  1.694e-05 3014.789   <2e-16 ***
## depth                 -1.890e-08  2.001e-07   -0.094    0.925    
## factor(area)2          8.785e-04  1.213e-05   72.410   <2e-16 ***
## factor(area)3          1.286e-03  1.578e-05   81.479   <2e-16 ***
## factor(area)4          2.208e-03  1.735e-05  127.278   <2e-16 ***
## factor(area)5          5.848e-04  1.215e-05   48.126   <2e-16 ***
## factor(area)6         -5.110e-03  1.624e-05 -314.678   <2e-16 ***
## factor(area)7         -1.993e-03  1.686e-05 -118.241   <2e-16 ***
## factor(survey)2013002  1.121e-02  1.043e-05 1074.952   <2e-16 ***
## factor(survey)2013005  1.196e-02  1.604e-05  745.567   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.08924    0.06015  -18.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.7884     0.1533   18.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.89428    0.04038  -22.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  13
##       Residual Deg. of Freedom:  661 
##                       at cycle:  20 
##  
## Global Deviance:     -4545.402 
##             AIC:     -4519.402 
##             SBC:     -4460.73 
## ******************************************************************

## 
## Family:  c("BCPE", "Box-Cox Power Exponential") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = Species_Hrs, family = "BCPE", data = spue.df) 
## 
## Mu Coefficients:
## [1]  0.06263
## Sigma Coefficients:
## [1]  -1.547
## Nu Coefficients:
## [1]  1.604
## Tau Coefficients:
## [1]  -0.445
## 
##  Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   670 
## Global Deviance:     -4253.19 
##             AIC:     -4245.19 
##             SBC:     -4227.13

The variability and relative difference in rank of species density per area and survey as observed in the table above was confirmed by the GAMLSS model (“Box-Cox Power Exponential” distribution, AIC: -4519, df: 4) where all surveys and areas were significant factors in the model. Depth, on the other hand, was not found to have a significant effect on the species density.

Functional diversity

We used the same approach as Stuart-Smith et al to exclude very few observations of a single species and excluded all species observations occuring in 2 or fewer stations (traps).

The ‘dbFD’ function only works when run in-line, not the whole code-chunk.

To avoid dbFD crashing ‘m’ was set in the ‘dbFD’ call (the number of PoCA axes kept as traits to do the FRic analysis). Several levels of ‘m’ were tested, and this could probably be tested further (to higher levels of ‘m’)

Results for n_spec > 2, m=10

Results based on the output from the dbFD analysis:

No area came out with the highest or lowest FD across all four metrics, but some patterns emerged:

  • Area 1 was highest in FEve and FDiv, but second lowest in FDis and RaoQ
  • Area 2 was consistently high in all metrics, and highest in FDis and RaoQ
  • Area 3 was above average in all metrics, except FDiv
  • Area 4 was lowest in RaoQ and FDis, and below average in all others
  • Area 5 was above average in all metrics
  • Area 6 was above average in FEve, but below average in all others, and had the lowest value in FDis
  • Area 7 was above average in RaoQ, had the lowest value in FEve, and below average in FDiv and FDis

Overall the highest Functional Diversity metrics were found in areas 1 and 2, while the lowest were found in areas 4, 6 and 7.

RaoQ is the most used metric (the one used by Stuart-Smith et al), and this is the metric we will use and discuss in the MS. It supports the analysis of species densities for Area 2, although Area 6 has the second lowest RaoQ which indicates that although the number of species is higher, their functional diversity is lower (or what we sampled was lower) compared to the other areas.

Stuart-Smith tests various RaoQs by taking out one trait at a time and this way ranking the traits according to R2.

## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 49 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4209541
Area FEve FDiv FDis RaoQ
A1 0.682 0.905 0.189 0.051
A2 0.637 0.849 0.263 0.079
A3 0.644 0.71 0.239 0.074
A4 0.614 0.745 0.2 0.048
A5 0.627 0.785 0.247 0.073
A6 0.665 0.65 0.179 0.048
A7 0.516 0.713 0.197 0.064

Rerunning the FD model excluding one trait at a time

Need to reduce ‘m’ to 9 to avoid crash:

Zero distance(s)Error in convhulln(tr.FRic, “FA”) : Received error code 2 from qhull. Qhull error: QH6114 qhull precision error: initial simplex is not convex. Distance=1.4e-16

Removing the MaxLenght trait from the model lead to the largest increase of R2 (by 0.047 compared to the model with all traits included), and higher than removing any of the other traits from the model. R2 was reduced when removing any of the other six traits from the model, indicating that all other traits should be kept in the model.

Evaluation of runs: - Use the model without MaxLenght (has the highest R2 of 0.468) - Removing Maxlength makes sence because the data included a few catches of sharks with very large MaxLengt (>2m) that skews the MaxLenght distribution. Excluding MaxLength creates an analysis around traits that are more common in distribution across species. - Also, the sharks are rowing predators, not so site attached (but that also applies to some pelagic species, especially the large ones, so the same argument applies for these - excluding large body size reduces the importance of these non-site attached species)

## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 37 PCoA axes (out of 46 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4684034 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 42 PCoA axes (out of 51 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3984771 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3663967 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3612455 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 48 PCoA axes (out of 57 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4029848 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3538643 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4132985
All MaxLength Trophic.Level Trophic.group Water.column Habitat Gregariousness diel
R2 0.421 0.468 0.398 0.366 0.361 0.403 0.354 0.413
All MaxLength Trophic.Level Trophic.group Water.column Habitat Gregariousness diel
A1 0.051 0.055 0.065 0.049 0.040 0.063 0.040 0.059
A2 0.079 0.102 0.103 0.059 0.059 0.093 0.070 0.089
A3 0.074 0.078 0.092 0.061 0.061 0.100 0.062 0.077
A4 0.048 0.063 0.063 0.050 0.031 0.046 0.039 0.057
A5 0.073 0.085 0.093 0.060 0.055 0.091 0.061 0.082
A6 0.048 0.050 0.058 0.036 0.042 0.064 0.046 0.051
A7 0.064 0.076 0.083 0.057 0.050 0.067 0.060 0.070

Results table

This table combines CPUE, traits pr gear type, species density and RaoQ into one results-table.

area TrapCPUE GillnetCPUE TrophicL_Traps TrophicL_Gillnets Greg_Traps Greg_Gillnets SpeciesDensity RaoQ
1 0.117 1.639 3.825 3.933 1.510 2.412 0.015 0.055
2 0.149 1.176 3.815 3.891 1.478 2.261 0.013 0.102
3 0.136 0.309 3.844 4.140 1.553 2.500 0.023 0.078
4 0.059 0.698 3.765 4.284 1.357 2.107 0.029 0.063
5 0.079 0.726 3.797 3.853 1.627 2.176 0.015 0.085
6 0.135 0.270 3.762 3.927 1.531 2.286 0.014 0.050
7 0.098 0.683 3.828 3.612 1.565 1.900 0.015 0.076